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conditions. Generally, for flood routing in rivers, the Saint-Venant equations will be used 

which can be solved by finite distinction method. Saint-Venant equations will be converted 
into nonlinear equations and will be solved using the Preissmann scheme in the finite difference 
method. Using the Newton Raphson method, the set of equations will be changed into linear 
equations and will be solved by the space method. Our aims are to the estimated zones of Saint- 
Venant equations for flood routing by using the finite difference method with over bank unsteady 
flow in an open channel. The effectiveness of this method to optimize the choice of finite 
difference method is more accurate than other methods having adequate space and time steps. 
Keywords: Flood routing; Saint-Venant equations; sprace method; over bank flow; hydraulic radius. 


|: this paper, we learn about the control of open channel water glide under the flood routing 
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Introduction 

Over the last two decades, the study of open channel water flow has become an active 
reseatch area due to its numerous advantages. It has been extensively investigated by many 
researchers. The hydraulic and hydrological problems involved in the computation of flood waves 
are based on Saint Venant (SV) equations, the pair of continuity and momentum equations. This 
study is important in an open channel under non-uniform or unsteady flow. It is worth mentioning 
that in this case SV equations cannot be solved analytically. Stoker [1] obtained the approximate 
solution of SV equations using the explicit finite difference method described in [2]. Ooi and Weyer 
[3] described the control design for an irrigation channel derived from physical data. Many 
researchers have obtained approximate solution of SV equations in a particular case of open 
channels ((4], [5], [7], [8], [9], [10], [11], [12], [13] and [14]). In these models, river waves can 
be categorize as gravity, diffusion, or kinematic waves, that relates to exceptional types of 
momentum equations. 

The flood routing is fully based on the unsteady flow (long wave or surges) and the water 
storage equations. At the point of the channel, a flood hydrograph will find out from the well- 
known hydrograph at various points upstream or downstream via the known channel characters 
with the characteristics of side inflow or outflow between these two points. Neuyen and Kawano 
[15] obtained simultaneous solutions for flood routing in the open channel network by using the 
Preissmann method for their proposed model. Kazezylmaz-Alhan et al. [16] discussed the reliability 
of the finite difference method for solving proposed diffusion and kinematic wave equations that 
describe the overland flow. Kohne et al. [17] presented the diffusion and kinematic wave method 
for runoff and surplus computation. Das [18] presented the Muskingum model to find the flood 
path and obtained its coefficients by employing the optimization method. Sulistyono and Wiryanto 
[19] investigated the flood routing by dynamic wave model in the trapezoidal channels. 

The main objective of this paper is to develop a quantitative method for identifying river 
wave types in the case of flood routing with the overbank flow. Different theoretical cases relate to 
the ratios between the central channel and flooded location concerning breadth and glide will be 
analyzed. Moreover, we define the estimated zones of SV equations for flood routing by using the 
Preissmann method with over bank unsteady flow in an open channel. 

Materials and Methods 

For flood routing issues like river waves and dynamic modeling of one-dimensional, the 
numerical answer of SV equations will be used. In the case of the flooded region in the river, let 
B, and By be the breath of the central channel and flooded region of the channel, respectively as 
given in Figure 1. Also let Ay and A, be the cross-sectional area of the central channel and flooded 
region of the channel, respectively. 


Figure 1. Open Channel having two breaths and cross-sectional areas for flow. 
SV equations consist of two equations; continuity equation and momentum equation. The 
continuity (mass) equation is given as 
dA, OQ 
= + — = q; 
Ot Ox 
which further implies that 
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Whereas the momentum equation is given as follows: 


aQ (aa) + ah 


ot 
lap aa) +gA2 6* + Sp) — gArSo + qv = 0, 
ot Ox 
Ay + Ap? +6A2 C+ Sp) — gAgSo + qV = 0. 
Thus, the momentum es becomes 
av av. ah qv 
aot t+e5, + 9(Sp — So) +3 = 0, (2) 


where is the flow high (m), is the velocity of flow (m/s), is the acceleration due to gravity 
(m/s*), is energy line slope, is the slope of the riverbed, is discharge (m*/s), is the horizontal 
distance (m), is time (s) and the q is a lateral inflow. 

To derive this system, a basic assumption will be made: that there is one-dimensional flow 
in both regions of the channel (central region and flooded region), and there is lateral inflow or 
outflow. 

Generally, the side channel will be rougher than the central channel. The velocity V in the 
central region is greater than as compared to the flooded region. In this case, the Manning formula 
can be applied separately to both regions in determining the velocity of both channels. In this case, 
the manning formula might be applied one by one to every section in determining the velocity of 
the section. After that, the discharge in the section will be computed. The total discharge will be 
equal to the sum of these discharges. Since velocity V in the central channel will be greater than the 
velocity in the flooded region so the component of discharge in the flooded area is small as 
compared to the discharge in the central channel; hence discharge Q can be approximated as 
follows. 


Q = AV, (3) 
where the cross-sectional area A, depends on x and t. As 
Ay = B, h. (4) 
Differentiating equation (3), we have 
2Q _ AY) _ Aaa, 4 av 
Ox ot _ Ox ‘ 1 ax 
_ By ov 
re a Ox = (Byh) Ox 
= Biba + (B,h) ae (5) 


Substituting equation (5) in equation (1), we get 


oh oh 
By + (BV) 5+ Bin) 


7 On as (v 72h 
2 Ot 1\" ax Ox a 
and so, 
Baty s ax thas (6) 


If 7 is the ratio of flooded region breath Bz and central channel breath Bo, that is, 7 = = 
1 


then we obtain 
Oh Oh , A1lOV 
ap ta a (7) 
ot Ox B10x B1 ; : 
Generally, in the momentum equation, the term will be calculated by the Manning formula. 


Since the velocity in the central region is greater than the velocity in the flooded region, thus the 
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term of the Manning formula in the flooded region is smaller as compared to the central region of 
the channel. As 


which implies that 


and hence 
4 
Se = (Vn)?R73 = (Vn)?R-™, 
where R is the hydraulic radius (m), n is the coefficient of roughness and mM is constant 


(m = =) and (h << Bj) is for large rivers. For central channel, the hydraulic radius is given by 
— Bh _ he 
= Byzh ye 0) 


Substituting equation (9) in equation (8), we get 
Spe nV-en™. (10) 
As Froude number is a dimensionless number which is used in hydrodynamics to specify that how 
a particular model works in relation to a real system, so 


V 
F=—, 


ign 


GhFo = Vv"; (11) 
Substituting equation (11) in (10) implies 
Sp =n? gF*ht™., 


that is, 


As 
OS 


—=n’gF7(1-—m)h™ le 
Ox 


Ox’ 

so, we have 

Oh _ nm OSf 
ax n2gF2(1-m) dx" 
Using (12) in (2), the momentum equation becomes 
OV OV ne OS ea qv . 
ae Oe age Ox + (Ss ota = , 
that is, 
av av hn™as V 
SHV tS + G(S-So)+Z=0, (13) 

where B = n?F2(1 —™m). Continuity equation (7) and the momentum equation (13) 
provide the general form of SV equations with flooded region and the side channels rougher than 
the central channel. In such case, SV equations depends on parameter that appears in the continuity 
equation. In the particular case 7 = 1, equivalent to flood routing with no overbank flow, was 
considered widely by Moussa and Bocquillon [12]. These two equations can be solved by using 
finite difference method. According to this method, Preissmann method is used to solve SV 
equations. 
Result and Discussion 

The Pressman model has a vast application in flood routing hydrograph in the open 
channel as given in Figure 2. 
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t 


z Az 2+ 1 


Figure 2. Preissmann method 


By using above model, fe time derivative is as follows 


0 
or lett = y [So va ae] os (1- w) js afi “| (14) 
Similarly, the space es is given as follow 
af |k+i _ [fit ia fits fi 
= | = o +(1—a) 4), (15) 
And the others terms are given as 
kt1, ¢k+1 kg ek 
f=a Je + (1-a) ese] (16) 
In the SV equations; for the terms ie and lie we will use equation (14), and for terms 
0 0 te) ‘ saa "= A . 
~, <r and - we will use equation (15), and for terms oe 2 = ve — and Sr, we will use equation 
1 

(16). The continuity equation (7), for p= - becomes 

hist — hia e hy? _ hi +I(a = + viet +(1-a) Via oF Vie x 

2At 2At Z 
A k+1 Ay k+1 
oth gy ths A] | cqy Buen * 60) 
Ax Dx . 
k 
(#) + Ve ag — vk ae vk 
+ (1 - @) 4+ a J. foo! ee += a) ee) 
Ax 
oo ae ay ey 
By i+1 By i By i+1 By i 
—4| (a) 5 + (1 - a) —— | = 0, 


which implies that 
A 
n[nieg + ak*t] — [hs + ak] +o? (visa! + VE) (nity — ni) 


At 
+Fa(l — a) (Visi! + VK) (his. — Ae) 
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At 
+ ie a(1—a)(VE, + VK) — Ak) 


tat a)? (V, ha + Vi) (hia — hi‘) 
k+1 At k+1 
Ste (af +(3) ) vittt — y+) 
k+1 k+1 
+ Salt a a((2) (a +(2) Jv i+1 — Vi) 


as By i 


rSaci-a)((*) + +(2) ) wm - vit) 
+a — a)? (2) + (2) ) (vk, — Vi‘) - qAta (¢) + =) ) 


i+1 i+1 


gar ae) (2). + (2) ) = 0. (17) 


Finally, a function T; for interval t will be obtained as follows. 
TM VEE) = an] + 


ee az (>) + (2) ykt1 _ ye += aD ykt1 _ prt 
Ax ( B, ve By ( i+1 ) o( i+1 l ) 


At A k+1 A k+1 At 
vies (GH) +), ) teas eer wry — a 
i+1 


At At 
+<-ab,(hky} - hi*t) + —aDe(V itt + VIF7)4D, — Dg 


— gAta (2) Es (2)”) =0. (18) 


In equation (18), D; to Dg ate ae which are given below: 


k 
—hiia - hj , 


D,=(1-@) o(@, (3) ) 


D3-(1— @)(Vih1 — Vi‘), 


At 
Daw py Pads 
Ds-a(1—- a)(Vi. +V*), 
Dee(= a)(hi — hr), 
At 
Dy = 5—DsDe, 
Ly 1)" 
ea) (G7 Gy) 
Similarly, the momentum equation (13), for p = = = becomes 
viet yk, yktt yk vt es VEL +VE 
| 2At — 2At +o 5 | 
yt — kt vE 
feos «aia 5) 
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1 hm k+1 + (h™ k+1 hm kL h™ k 


2 (Se = (5), -) ae (Se = (5), } 


x 
Ax Ax 


_ (cai : (s;) . ace (Se (s;), } 


V k+1 V k+1 V k V k 
a Gis at) wd a" 7), | _ 


—gSo +g 


which implies that 


A A 1 
ong WES + VE] — [iia + VF] + atu + ve) (vi — vi) 
1 1 
A a(t — a) (Via + Vi) (Vist — Vi) + G (1—@)?(Vii + Vi)(Vih — Ui) 


an 
gP 


¥ +1 1 
«(SY = AE) + aa = (ernest + EM) (SIE, - GH!) 


: + +1 1 
+ aa (Mba + OD) (SIE GI") +3, a0 


(HE + HE) (Sh = GE) - 950 + dea (SEY + (SIE) 
: : gAx V k+1 V k+1 
+ Ax(1— a) ((S;),,, + (Sy); ) + a (2) + (x) 


i+1 i 
k k 
qAx V V 
*(1-a)((4) +(4) )=0. 
+7 COG). Fe), (19) 
Consequently, a function U; for interval t will be obtained as follows. 


Ax 
UCR VE hist View) = Tag [Vier + Ver") + Do 


1 
45 a(1— a) (Vita + VE) (Vie Vi) + Tea? ((h™ Diet + (HE) 


2 
i . (visit) r (vi) + Dio (Vit? — Vi") + Dis (Vist? + Vi**) + Dip 


1 + +1 
+ ar(nm sist + HE) (Sieg — (SG), |) + aia (HORE + DE) 


1 
+75 aD (SS) = (s;),"’) + Dig t+ Dig t Axa (Ss) 4 (s;),"’) +Diy 
qAx V k+1 V k+1 _ 
+2 a (oe + =), ) + Dig = 0, (20) 


where Dg to Dyg are given as: 
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‘ 
Di ri a(1—a@)(Vii, —V*), 


Diy = 20? | UE) - (vi)"], 


1 
D43 = ae = a) (Com 7 (s,);), 


Dy = (1 - a)((h™ 4 + (ae); 
Dis = >; Dy4D43, 
Di6 = —2AxSo, 


Dy = Ax(1—@) (Con - (5;);) 


ee am ma-a((Z 4). +()) 


As in favor of each interval, - two vector functions and ate the functions of four variables, 
namely. In the favor of each interval, we define two equations, and similarly, for intervals equation 
will arise. For each interval, we have a node. And for each node, there will be two unknowns (flow 
rate and depth), thus if there are unknowns, then two equations are formed. Thus, we will obtain 
two more equations from up and downstream conditions [20]. 

Finally, a system of nonlinear equations will be generated; to change the nonlinear equations 
to linear, we will use the Newton-Raphson method. According to this method, the derivative of 
functions is given as follows 


k+1 
[cankt] - ai (Fe), 


oT; A k+1 k+1 
SuK+1~ 77 avid —V; 
dhe* | i+1 i | [ep 
dB 
At [ek] - Ar yet (dy At oryktd _ kt 
ie 83 |] a it — 0] 
Ax [ee] Ax 
1 
GH), 
——— aD, + gqAxa : 
A 5 (pak 
that is 
kt1 aBy k+1 
= = n+ elvis =v _ ae Ga), |, At aD 
a. K+1 i 5) 
Oh; [Bk] are 3 
_ agen(eay 7 . vate sae 
[(By ey eres a [vist - Vit] —~ Ax aDs + | (F), . (21) 


Also 


Dec 2021 | Vol 3|Issue 5 Page | 103 


OPEN ACCESS S % 2 = 
3 International Journal a ee in Science & Technology 


dA dB,\*** 
k+1 1 k+1 1 
oT; At k+1 _ yk+1 Bviva Cia i, 7 Aviv Cae 
peer = +50 [Man — Vi je 
oni [eds 
k+1 k+1 
dB 
k+1 1 
At (BST (Sa Th), — Aditi (TE). At ace at 
te Oa Na toe ar | 
. [Bis . 
dB, k+1 
At (Fe), 
+— aD, + qAxa——,, 
Ax [@pey” 
that is, 
dB k+1 
A,)Kt1 (761 
ie ~ 2fyk+1 _ yet OO a a 
k+1 n+ a Vi+4 i | k+112 
one [vi 
dBi\ 
k+1 il 
At (Aina At 
+a as 1 _ GON Cai oa) [vit — vk] += ads 
7 [Bi = 
_ qhta _ (=) k+1 (22) 
[(By ykeay? Jeet 
fos we te 
aT; AG: y (2) mi At At 
ey ee BB pad ——aD +e’ hkt1 — nkt+1 
one Ax By eg By Ax 2 [ i+1 i | 
+= aDe, (23) 
at; At 9 ffa\**t . fA ) At 2 {hktt — pk+1] 4 ot 
ane Ax (Ce J ee ax *P2 + [nist — hi lear aDe, (24) 
k+1 
aU; _ a ia Rel “aye _ © +4 k+1 a 55), 
ohn = Fe a (Sp), mca di - aoe ty Ch); ) ~anktt 
k+1 k+1/dAz\K+1 
1\k+1 aDy4 a(sp),* a(sr); vin (Ge ), 1 
+aD;3(m(h™")p*") + one -+Axa one ~Tankay (25) 
As Sp =n 2gF7h'-™, we oo 
a(S), 2472 —m)k+1 
ohkt1 =n gF (1 —m)(h )i , 
U 
that is, 
k+1 
a(Sf). 
ya = gBA™eE, 
U 


ve B = n’F*(1 —™m), and the above equation (25) becomes 


OU; k+1 k+1 = - 
aneH = gB La Sica = (S;), mca ile = a? ((n™ yt + (am) )(A a 
k+1 k+1 
D pm-1yk+1 D..(h-™) +1 + A p_rmyk+1 _ Vina dA 26 
+aD;3(m( )Er*) + aDyg(h-™)F** + Ax agB(h-™); aCe (26) 
Also 
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Ou = a? ((s ee = (Ss *) many eo a (cami ait ane a(S), 
ane gh Ti ne nape Uae ane 
k+1 
k+1 k+1  y7k+1 (GA2 
m-1\k+1 aD, a(S,), a(Sr), Vita ( dh Di 
+aD,3(m(h Vivi + — er ae aie ee ae 
gB Oh Ohisy (Ca 
that is, 
aU; _ a ile ies m-1)k+1 4 ,2(¢p,myk+1 myk+1)(p—myk+1 
Spe ~ gp (Shinn ~ (Sr) mE + a2 (HME + GEO MEE 
= _ 2 vit da,\Kt1 
+aD13(m(h™ 1) ) + aD (h-™ i + Ax agB(h-™ it — aac ‘ (27) 
27i+1 
In the similar way, the partial derivatives of remaining equations are 
OU; Ax 2a? qAxa 
aver = ar ae — Pro + Dia + Fat (28) 
and 
OU; _ AX | 2a? 4d qAxa 
avEHT = Gat are Vi + Dip + Dy + Taper (29) 


As a result, in support of interval, equations can be produced. Hence, it is ample to get 
the other two equations from the upstream and downstream boundary conditions, so for all nodes, 
two unknowns and values will be obtained. The upstream borderline (node No. 1) may 
additionally be the inflow hydrograph, where is the discharge of influx hydrology, and the 
downstream boundary circumstance (node No.) may be the glide phase, the place is the overflow 
height, is the overflow coefficient factor, and is the weight. 

Thus, the downstream and upstream conditions are given a 
(hy, Ve) @ © (hat Viet) 


(hy, V,) @&——_— ® (hy, V2) 


As 
To (hy, Vy) = Q1 — AV, (30) 
Ty41 Ast Vert) = AnstView1 — CW(Aey1 — HY. (31) 


The initial condition is the steady flow before the flood, so the partial derivative of equation 
(30) and (31) are 
ATo(hiVs) OQ: _OAVs __, OAs 


pees ae ea 32 

Oh, dh, Oh, Ona oe) 
OTo(h1,V,) OQ, OAV, OV, 
eS ==), — = =/ 33 

AA av, «av, "AV, - oy 

OT 41 Asa Vier) _ OAx+1View1 IC C.W (hear — Hy?) 
any ah, ah, 
aA 
Ves i — (1.5)CW (hes — HY, (34) 
and 


OT 41 Nest Viet) — OAx+1View1  ICC.W (hear = Hy") 
OVi+4 OVi+4 OVk+4 
= Arti (35) 
To solve non-linear equations from (32) to (35), we will use the Newton-Raphson 
algorithm. 
The Newton-Raphson method is one of the most familiar iterative schemes used to solve 
non-linear equations. First, we will write given equations in vector form: 


T,(%1, X2, ag oa) = 0, 
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where i = 1,2,3,...,2N,X = (1, Xp, ..., X2y)denotes the vector of unknown variables. 
By Taylor series expansion, we have 


“aT; 
T(x + 6x) = T(x) + = =" 5x; + O(5x)? 
iz. 


fori = 1,2,3,...,2N. In the above equation, first partial derivatives form a Jacobian 
matrix. Ae we take 


c= yz 


then 
T(x + 5x) = T;(x) + J 6x; + O(5x)?. 
If we neglect higher order terms and set the left hand-side equal to zero, then we obtain a 
set of linear equations given by 
J 6x; = —T;(x). 

The above system of linear equations in matrix form can be further solved by using 
Gaussian elimination method or LU decomposition method for the unknown value of x, and hence 
the approximate solution is obtained by 

Xnew = Xoia + OX. 
The iteration procedure will be carried out until a reset convergence stage is achieved. Now 
in our case, we write our above-mentioned equations in the same way as above 
Ty (ht **, V**) = 0. 
By Taylor series expansion, we can write as follows 
OT) 
To (htt? Vt? Ant, AVI*1) = Ty (hy Vi.) + aaett Aa ae 


Pot The Vi an AVE) =0 and let T)(hy,V,) = = tg, so the above equation 
will become 


AVi+t, 


OTo OTo 
to + ppp Ant + apna Avi'** = 0, 
that is, 
OT» OTo 
ane ant + Sper AVE = ty. (36) 
Also, 


T; (Cee Ve} = 0, 
By the Taylor series expansion, we have 


T, eV an Ae = T, Aye + Ahi? + 
l l U U l L l U 


ane ae 
ae 


ayir 


1 
ayaa Ah}* 5 avr 


i+1 
Put T,(hp**, Vit" Anpt* A cae = 0 and 7;(h?'*1,V;"*") = t;, we get that 
aT; aT; 
+1 +1 +1 +1 
anett aa any + we AVP et aaa i eres avnst Av," =). (37) 


i+1 
Similarly, we obtain the following equations 


au; au; 
1 1 1 
ane pal Ahj** + we AVtT + 4 Apts Ahj** + ae AV?+* = —uj. (38) 


i+1 i+1 
OTn+1 1 OTy+1 1 
anntt i+ avert AVNSi = twat (39) 


+ 
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In the above equation (36) to (39), to, tj, Uj and ty41 be the values of Tp, Tj, Uj and Ty 41, 
respectively. Finally, we get a system of 2N + 2 equations and 2N + 2 unknowns, where 
unknowns ate h and V. The system of equations are transformed into the matrix form is as follows 


memo... 1... 0 Ahy a 
aT, OT, OT OT. AV, a. 
Ohi OVi Ohg OVa2 : Aho :! 
BU, BU, AU, OU; AV, = 
Oh, OV, Ohg OVo2 * 1 : Al —ug 
0 Q 222 272 aT2 aTr 13 

Oh2 OV2 Oh3 OV3"°— * ; Av3 |= 
0 0 OU2 OU,g OU2 QU2 F 

Oh2 OV2 Ohg OV3° 
0 OTN b1 OTN hl Ahn+1 —UN 

“Ohn+1 OVN +1 AVn+1 —tn+1 


The mattix elimination programs, such as Gaussian elimination or LU decomposition method 
will be used to solve the above given system. However, Fread [9] used the Sparse characteristics of 
the Jacobian coefficient matrix with maximum value are considered for consecutive elements and 
developed an effective algorithm to solve such kind of Ribbon matrix problem. No matter which 
solution is used to solve the matrix, it is ultimately the case. The step is to calculate the correction 
value of the unknown of the next iteration. Thus, we obtain the following 

Qi*" = QF + 6Qi, 
Att = hf + Shy, 
Numerical Experiment. 

The above method requires main and boundary conditions. The main condition is a stable 
flow before the flow change, where the upstream condition can be flooded hydraulic pressure, and 
the downstream boundary condition can be a displacement relationship. 

We consider a rectangular open channel with a length of 3 km, a breath of 3 m, having bed 
slope 0.0005 and a Manning roughness coefficient at the beginning of the channel is 0.25. When 
the flood entered with the specified water level (in flow water level) with a spillway at its end having 
height of 1m, width of 15 m, and coefficient factor 1.6. Then, the water level curve and flow water 
level on each point will be drawn at any time by the above method. For calculations, MATLAB 
software has been used. 

In Figures 3 and 4, the flow of water graph at 500m and 1500 m in intervals from the 
flood beginning point has been drawn with Preissmann scheme. According to these Figures, we 
observed that through the increase of distance from the starting point of flood, difference between 
results increases. 


xz=500m ams FD | 


PREISSMANN METHOD 


120 180 


Time (3s) 


Figure 3. Hydrographs computed at the 500 m distance from the flood start 
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Figure 4. Hydrographs computed at the 1500 m distance from the flood start 


Conclusion. 


We study the control of open channel water flow under flood routing conditions by 


employing the Saint Venant equations. The Saint Venant equations converted the structure of flow 
in the given channel into nonlinear partial differential equations which can be solve numerically by 
using numerical methods such as finite difference method and Newton Raphson method. By using 
this method, we develop iteration scheme to estimate the flow and height of water in the given 
channel. We also study the estimated zone of Saint-Venant equations for flood routing with over 
bank unsteady flow in the open channel. Some numerical experiments are also presented. 
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